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THE APPLICATION OF THE ROUTH APPROXIMATION METHOD TO TURBOFAN ENGINE MODELS 
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National Aeronautics and Space Administration 
Levis Research Center 
Cleveland • Ohio A4135 


Tl;e Routh approximation technique is applied in 
the frequency domain to a 16th order state vari- 
able turbofan engine model. The results obtained 
motivate the extension of the frequency domain 
formulation of the Routh method to the time domain 
to handle the state variable formulation directly. 
The time domain formulation is derived and a new 
characterization, which specifies all possible 
Routh similarity transformations, is given. This 
characterization is computed by the solution of 
two eigenvalue -eigenvector problems. The applica- 
tion of the time domain Routh technique to the 
state variable engine model is described and some 
results are given. Additional computational prob- 
lems are discussed including an optimization pro- 
cedure which can improve the approximation accu- 
racy by taking advantage of the transformation 
characterization. 

INTRODUCTION 

In control system design the system to be 
controlled is often represented by a complex math- 
ematical model. Practically, this complex model 
may be difficult to use for design purposes. 
Additionally, the resultant control design may be 
too complex to implement. To eliminate these 
problems model reduction methods are often em- 
ployed before designing the control system to re- 
duce the complexity of the original model while 
maintaining the important characteristics. 

One philosophy adopted in reduction methods 
is dominant mode approximation (ref. 1). in this 
paper one such dominant mode approximation tech- 
nique, the Routh approximation technique (refs. 2 
and 3), is applied in both the frequency and time 
domains to the simplification of a state variable 
turbo fan engine model. 

Hutton and Rabins (ref. 4) have shown the 
frequency domain formulation of the Routh approxi- 
mation technique to be an excellent reduction 
method for many complex mechanical systems. This 
formulation is applicable to a turbofan engine 
model as well and is therefore a strong motivation 
for the work of this paper. 

The objectives considered include the appli- 
cation of the Routh approximation method in the 
frequency domain to the FIDO turbo fan engine. 

This application demonstrei^s that adequate, re- 
duced complexity, engine models can be determined 


via the Routh technique. Next, the frequency do- 
main Routh formulation is extended to the time do- 
main. This extension allows the usual engine time 
domain formulation to be handled directly. The 
time domain Routh approximation formulation also 
allows additional flexibility in the choice of the 
approximate system. Using the additional flexi- 
bility of the time domain formulation, optimized 
approximations may be readily generated. The time 
domain formulation is applied to the engine model 
and some results given. Approximant optimization 
results and the application of the Routh method 
for deriving fixed dimension realizations of the 
engine example are also discussed. 

ROUTH APPROXIMATION 

Modern jet engines can be represented as mul- 
tiple input -multiple output (MIMO) dynamic systems 
using linear, constant coefficient differential 
equations to model small variations at various op- 
erating point conditions. A transfer function 
representation of these equations can be used as 
input information in a classical or multivariable 
frequency response control system design. The 
complexity of these transfer functions can cloud 
the design process with unnecessary detail or dif- 
ficult computational problems. Additionally, con- 
trol system complexity may be adversely Influenced 
by unnecessary model complexity. Thus, there is 
a strong motivation to reduce the complexity of 
system and control models. 

Frequency Domain 

The Routh approximation method (refs. 2 and 3) 
is a dominant mode reduction technique. In the 
frequency domain the method incorporates parameters 
obtained from a Routh stability table analysis in 
the reduction process. These parameters, called 
alpha and beta parameters, are obtained from the 
coefficients of the MIMO transfer function 
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where the a^ are scalars and the are f x m 

matrices. Algorithmically, the n alpha param- 
eters are computed from the denominator of (1) as 
(ref. 4) 
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where 


■ 0 for i -*• J > n 

The coefficients a^ can be arranged as entries in 

a Routh table (ref. 2). This tabular construction 
suggests the name Routh approximation. 


The beta parameters are computed from both the 
numerator and the denominator of the given transfer 
function. Consider the single input-single output 
(SISO) transfer function from input q to output 
p implicitly given in (1), Hpq(s). The beta 
parameters are found for this SISO transfer func- 
tion as 
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where 
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for 1 + j > k 


Note that for different output-input pairs differ- 
ent approximant orders, k, could be chosen. The 
k-th order approximant for the pq-th output-input 
pair is written as 
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bj - 0 for i + j > n 

and bp is the qp-th element of the matrix 
Bp, r - 1, n. In this way beta parameters for the 
i X m possible SISO transfer functions of (1) can 
be found. Again considering the SISO transfer 
function, Hpq(s), the computation of the k-th 
approximant from the first k alpha and beta 
parameters is given recursively as 


(8) 
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Properties 

The Routh approximation method exhibits sev- 
eral very useful properties. Each is important in 
the reduction problem and is described briefly. 

Stability . If the origins, transfer function 
is asymptotically stable, the Routh approximant, 
of any order, will be asymptotically stable. 

Pole-zero locations . The poles and zeros of 
the approximants approach rhe poles and zeros of 
the original function as the order of the approxi- 
mation is increased. 


Impulse response energy . If kM is the im- 
pulse response of Hpq (s) , and if the "impulse re- 
sponse energy" is defined as 
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Alto, if ^^(t) It the impulse response of the k-th 
order tpproxlmant of Hpq(s) , then 
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Since ^q(^) 1^ assumed asymptotically stable, 

> 0 for i ■ 1 n and 

0<E, <E, <£-<...<£ -E (13) 

12 3 n 

The ratio gives an indication of the approx- 

imation accuracy in terms of the percentage of 
total energy accounted for by the k-th approxima- 
tion. 

Derivatives . The k-th Routh approximant sat- 
isfies the following derivative condition. 


The system is assumed stable with distinct 

eigenvalues. For the jet engine systems studied 
here this is a realistic assumption. Note that 
the reciprocal transformation of £2 is 
This transformation will always exist since A is 
assumed stable. The transformation is required to 
preserve the dominant or low frequency information 
of Zi during the reduction process (ref. 3). 

Now consider a nonunique similarity transformation 
of £2f called the Routh transformation, such that 

ATp - Tj^R (18) 
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where R is the Routh stability matrix. That is 
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(19) 


Initial and final values . The initial und 
final values of the step response trajectories for 
the original system and its k-th order Routh 
approximant are the same. Note here that this is 
accomplished even though the approximant is a 
strictly proper transfer function for all p and 

q- 


Time Domain 

The Routh approximation method is being ex- 
tended to the time domain to enable Routh approxi- 
mants to be directly determined from a state-space 
model formulation. The alternative would be to 
determine the p ^ q transfer functions for the 
MIMO system and then reduce each function individ- 
ually. For a typical engine model with p - 16 
and q * 3 this alternative would involve the re- 
duction of eighty transfer functions at each oper- 
ating point. Computationally, this would be unde- 
sirable. Hutton (ref. 3) outlines a method for 
directly determining the alpha and beta parameters 
from the state-space formulation. However, compu^ 
rationally, this method failed for the 16-th order 
engine example of this paper. Thus an alternate 
computational procedure was devised that incorpor- 
ates eigenvalue-eigenvector solutions which are 
well known and for which excellent computational 
solutions exist (ref. 8). 

Given the state-space representation, 

• A, A 
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and 



- diag[a^) i - 1, . . . ,n (21) 

Also, consider the modal transformation of Z 2 
where 


AT D - T D A (22) 

m m m m 

and 

A - diag[X^l i ■ 1, . . . »n (23) 

H<tre the are eigenvalues of A and the i-th 

column of Tqj is the corresponding eigenvector. 
Note that Tj^ will be unique with respect to a 
scaling convention and a specific order of the 
Xj^'s. The scaling convention is represented by 
an arbitrary, full rank, diagonal matrix with 
diagonal elements d^ that may be complex. 

Since R is similar to A 


where xeR^, ucR°^, and ycR^, define a reciprocal 
system, £2 such that 
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Again la a laodal matrix for R and A, and it 

ia uniqua with reapact to a acaling convention rep- 
reaented by Dz* From (18) » (22) • (2A) • and the 
definition 

D - (25) 

n _ 

all poaaible Routh tranaformationa can be charac- O 
terized as 

T_ - T (26) 

l\ m '■r' 

Recall that a fixed ordering of eigenvalue*: ia 
given. Permutationa of thia ordering can be repre- 
aented by a permutation matrix, Ep, where 


Examination of (32) ahowa that the reduced 
order ayatem matrix, ia aimply a aim- 

ilarity trana format ion of the truncated Routh ata- 
bility matrix. A truncation of an n-th order 
Routh matrix (R) yielda another Routh matrix of 
k-th order (R^i)* Since the original ayatem ia 
atable, the alpha parametera of R, a^, i ■ 1, n 
are positive. Therefore, ia alao stable. 

Add iticr illy, the reduction process implied in 
(3?) exhibits the same pole properties as the fre- 
quency response approach. 

* -' Hutton (ref. 3) has shown that if the system 
given by (16) is a single input-multiple output 
system and Tg is selected such that 


Trp “ Vr 

and Ep is the identity matrix with appropriate 
columns interchanged. Although this permutation 
does not affect the input-output transfer relation- 
ship of the original system, it will affect the 
input-output characteristics of the approximate 
system. Undoubtedly, one of the n! possible 
state permutations will yield a better approxima- 
tion in some sense than another, but this problem 
ia not considered here. 

From (18) a Routh canonical system, £3, can be 
written as 


^ • «*R + 
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then the elements of the rows of H, h^, become 
the beta parameters for the output considered 

\ “ ^®li»®2i ®ni^ 

and the properties of the approximation described 
previously hold. It can be shown that ii* D is 
selected as 
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The Routh approximation procedure starts by assum- 
ing that 
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porating the ass\xnption of (30) becomes or 
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where 



and the subscripts 


indicate conformable partitioning of the appropri- 
ate matrices. The desired approximation then is 
the reciprocal transformation of called 


then (33) is satisfied. 

COMPUTATIONAL ASPECTS OF ROUTH APPROXIMATION 

Computational aspects of both the frequency 
and time domain formulations are discussed in this 
section. 


Frequency Domain 

In the frequency domain the Routh approxima- 
tion problem is solved by the calculation of the 
alpha and beta parameters for a given transfer 
function. Straightforward programning of equa- 
tions (3), (5), and (6) gives these parameters and 
the Routh approximants. This Routh table proce- 
dure is computationally efiicient and accurate. 

For example, it has been employed directly to re- 
duce a 43rd order line dynamics problem (ref. 9). 
Also, the frequency domain engine results were ob- 
tained using this procedure with no computational 
difficulties. More complete engine results using 
the frequency domain formulation are given in 
ref. 9. 


Time Domain 

Computationally, given the state-space formu- 


latlon of (15), cIm domain Routh approximation re- 
quirca (a) the alpha coefficients and (b) the Routh 
transformation matrix. The alpha coefficients can 
be found readily from (3) if the system character- 
istic equation coefficients can be found. A pro- 
gram that incorporates the method of Danielevski 
(ref. 6) t? find the system^characteristic equa- 
tion of the system matrix, A, was used in this re- 
port. 


A method for the computation of the Routh 
transformation matrix of (18) has been proposed by 
Datta (ref. 7). This method was programmed and 
applied to the 16th order jet engine model. Nu- 
merical results were unacceptable due to large com- 
putational errors and another technique was sought. 
Since accurate and efficient eigenvalue-eigenvector 
techniques are well known (for the case of distinct 
eigenv^ues), the computation of the Routh trans- 
formation matrix was reformulated as the solution 
of two eigenvalue-eigenvector problems. These two 
problems are represented by (22) and (24) and 
solved using the computet - t, methods of refer- 
ence 5. Since A and R are similar, they have 
the same eigenvalues. The Routh matrix R is also 
known from (19) since the alpha parameters have 
been found. 

The transfc rmation matrices of (22) and (24) 
may have elements in the complex field. Thus from 
(2) to (6), T{( can in general be complex. However, 
a real Tg is desired to yield a physically real- 
izable approximation and to facilitate computer 
computations. The matrix, Tg, can be constrained 
to be real by the proper selection of the elements 
dj^ of D. To determine these elements consider 
the eigenvector solution of (24) using real number 
computer operations. Given A and its eigenvalues, 
the modal matrix of A can be written as 

T -TV (35) 

nun mm 

where is block diagonal and Tj^q, the modified 

modal matrix of A, is a matrix of real numbers. 

If the eigenvalues of A were all real, then 
Vjj - I. HoweveiT, if and X^+l are a complex 

conjugate pair, then the i-th block of is 

defined as 
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for (24) 
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Now 

from 

(26), (33), and (35) and the 

definition 




(38) 

Tr 

can 1 

be written 




T^ - T D„T’^ 
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Now if Tg is to be real, then Dg, a block diag- 
onal matrix, must be real. This can be assured by 
selecting the elements of D, d^^, such that: 


(1) d^ is real if X^ is real 

(2) d^, are a complex conjugate pair if 

Xf^Xi^l ^ complex conjugate pair. 

It is remarked that the reciprocal transfor- 
mation of (17) does not represent a significant 
increase in computations since the modal transfor- 
mation is the same for A and A and the eigen- 
values are simply reciprocals. 

The matrix inversion of (26) can be elimi- 
nated by noting that 

r ■ ir^T (40) 

where 

T - dUgl(-l)^’^] (41) 

From (19) then 

AjjTR - R^AjjT (42) 

Now from (24) , (26) , and (40) 

Tp ■ T„D tJtA^ (43) 

K m m z K 

APPLICATION TO ENGINE MODEL 

The frequency and time domain Routh approxi- 
mation formulations were applied to a 16th order 
linear state variable representation of engine dy- 
namics for small deviations about an operation 
point. The actual engine has five inputs and six- 
teen outputs when the states are taken as outputs. 
However, for the purposes of this paper only two 
inputs and two outputs were considered. The 
inputs are 

Uj^ - w^ - Engine fuel flow 
u^ * A^ ■ Engine exhaust nozzle area 
and the outputs are 

y^^ - - Engine compressor rotational speed 

^2 * ^T * turbine inlet temperature 

These variables were selected because of their im- 
portance in engine control applications and be- 
cause of the wide range of dynamic content of N^ 
and Tj. In an engine the speed, N.» is one of 
the slowest responding variables while the temper- 
ature, Tx, is one of the fastest. The input data 
for this engine are given in Table V. 

Frequency Domain Application 

As an initial experiment the frequency domain 
formulation was applied to two 16th order SISO 
transfer functions. These transfer functions rep- 
resent dynamics of the turbofan engine from the 
input, W£, to the outputs, N^. and Tj. The coef- 
ficients for these transfer function were calcu- 
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lAt»d I'Jiing the progreme of reference 9 and are 
given in Table I. Table II shows the impulse re- 
sponse ratios for Routh approximations of increas- 
ing order for the two transfer functions. These 
ratios were calculated from the alpha and beta co- 
efficients as outlined in (12). From these ratios 
an estimate of an acceptable approximation order 
can be obtained. 

The approximant order was selected by first 
choosing a minimum acceptable level of accuracy 
defined by the impulse response energy ratio. The 
level chosen corresponds to a ratio of 0.81. This 
level of accuracy was assumed adequate for the 
purpose of this paper. Next the order that most 
nearly corresponds to the selected ratio was chosen 
as the order of the Routh approximant. For the 
energy latios giv*5n in Table II the Routh approxi- 
mant orders as founo by the criterion described 
sbove are 

\ - 3 (y^ - N^) 

- 9 (y2 - T^) 

The Routh approximants for and T.^ were 

calculated using the approximant orders and 

k 2 and the algorithm given in (6) to (9). The re- 
sults are summarized in Table III. Comparisons of 
pole locations and step responses for given and 
approximate transfer functions were made to evalu- 
ate the adequacy of the approximations. Table IV 
gives the pole locations and the step response com- 
parisons are given in Figure 1. Note that the pole 
comparison shows a good correlation between actual 
and approximate system poles. Likewise the step 
responses of Figure 1 indicate excellent agreement 
between actual and approximate model representa- 
tions. Rased on these results it is concluded that 
the Routh approximation technique is a viable ap- 
proach to the reduction of complexity In frequency 
domain models of Jet engine dynamics. Also, the 
accuracy level selected in this initial erample 
may be too stringent based on the step response 
comparison. Further results are summarized in ref- 
erence 9. 

Time Domain Application 

The time domain calculation procedures out- 
lined above were applied to the engine model assum- 
ing that a 5th order approximation for each input 
was adequate. A comparison of alpha parameters 
calculated by the frequency and time domain ap- 
proaches summarized in Table VI shows the accuracy 
of the time domain approach formulated in this 
paper. Step response trajectories for the original 
and approximate systems are shown in Figure 2 and 
indicate, in the case of T>p, the tradeoff between 
accuracy and complexity. 

DISCUSSION 

Note that for a five input system approximated 
by 5th-order models, a total reduced system reali- 
zation may require a state vector dimension of 23. 
In fact the realization dimension will probably 
vary for different numbers of inputs and orders of 


approximation. Often, however, the allowable 
dimension of the total reduced approximation is 
fixed by some engineering or economic constraint. 

Constructing a fixed dimension realization 
oT the total system from the Routh approximation 
can be a very difficult computational problem 
especially if the numbers of inputs and outputs 
are large. Thus the technique as outlined does 
not directly handle the fixed dimension problem. 
However, if the same Tg matrix is used in the 
time domain reduction process for each input, the 
total approximation can be realized by a system of 
k-th order where k is the number of alpha param- 
eters retained and the fixed dimension approxima- 
tion problem is solved. The selection of the Tg 
matrix which will give the best k-th order approx- 
imation is therefore an interesting problem. 

In this legard the somewhat arbitrary selec- 
tion of D can be used to good advantage. Many 
different transformations can be found very 
quickly once the original two eigenvector problems 
have been solved. Indeed, the selection process 
could be automated by performing an optimization 
of some function of error between the original and 
total approximant system over the n parameter 
space of D. Such an optimization scheme was 
tried on the example given in this paper. Two 
different error functions were used for compari- 
son. The first was a weighted sum of the differ- 
ences in system and approximant step response 
energy. The second was a weighted sum of squares 
of the differences in system and approximant 
steady-state values. Significant minimization of 
each error function was easily achieved using a 
conjugate direction optimization scheme in rela- 
tively small amounts of computer time. Thus, the 
general optimization procedure appears to be a 
good approach to improve the accuracy of approxi- 
mants for linear systems while maintaining a fixed 
order of realization. However, the time domain 
Routh approximation procedure when constrained to 
yield a fixed dimension realization in the mul- 
tiple input case does not exhibit the final value 
property of the single input case. For the engine 
example posed, the significant improvement in the 
approximation gained via function optimization was 
overshadowed by these final value errors for the 
multi-input case. 

Based upon this observation the time domain 
formulation was modified to insure that the final 
value property would be met by a fixed dimension, 
multi-input Routh approximant. In the original 
formulation the difficulty with fin^^ values in 
the multi-input, fixed dimension pro. cm can be 
traced to the original assumption of the reduction 
process (30). If the initial assumption were 
changed to 

Xp - 0 (44) 

then the approximation would force the final value 
property. Unfortunately, the assumption of (44) 
when applied to the original system may yield an 
unstable approximant for a given, stable system. 
This was the case for tne engine example and, 
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thus, ths Bodificstlon of tUs Routh procedure was 
rejected. 


CONCLUSIONS 

In this paper the Routh spproximstion process 
vss refonmilsted in the time domsin. A new chsrsc- 
terisstion of the nonunique Routh similarity trans* 
formation was derived which describes all possible 
Routh transformations. This characterization casts 
the computation of the Routh transformation into 
two eigenvector-eigenvalue problems which are 
easily solved. The application of the time domain 
formulation to a 16th order state variable descrip- 
tion of a turbofan engine was described a.id results 
given. These results indicate that the time domain 
Routh approximation technique can be a valuable 
uecnnique for reducing engine model complexity when 
dealing with the isodel on a single input basis. An 
optimization procedure was discussed that can sig- 
nificantly improve the approximation in a computa- 
tionally efficient manner by taking advantage of 
the new time domain characterization. 
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TABU 1. COEFFICIENTS FOR TWO 16TH ORDER 
ENGINE TRANSFER FUNCTIONS WITH 
FUEL FLOW INPUT (w^) 


TABU 111. ROUTH APPROXIHANT COEFFICIENTS 
FOR TWO TRANSFER FUNCTIONS 


Ng (P ■ 1. q ■ t) 


T^. (P ■ 2, q 


1 ) 
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»c<P ■ !• 
q ■ 1) 
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1.0000 


1.0000 
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2.6132 0.43978 


4.7092 

3.6223 





2 

2.5702 .61711 


1.0446E03 

9.3469E01 

0 

1.0000 




.8343 .25676 


1.3858E04 

1.3400E03 

1 

1.0638E03 

0.1140 

5.7270 



1.1722E05 

1.1906E04 

2 

3.7805E06 

4.2916E01 

5.8161E03 



6.3304E05 

6.6700E04 

3 

6.6915E08 

-4.2414E03 

1.8835E06 

i ■ 6 



2.1019E06 

2.2877E05 

4 

7.0216E10 

-2.2082E06 

2.9209E08 

7 



3.9634E06 

4.4718E05 


4.7772U2 

-8.2163E07 

2.6079E10 

8 



3.6848E06 

4.3229E05 

6 

2.2150E14 

2.2542E10 

1.4714E12 

9 



1.1961E06 

1.4512E05 

7 

7.1949E15 

2.8702E12 

5.5192E13 






8 

1.6584E17 

1.5789E14 

1.4142E15 






9 

2.715m8 

5.0230E15 

2.5027E16 






10 

3.1255E19 

1.0028E17 

3.0507E17 

TABU IV. COMPARISON 

OF 

POLE LOCATIONS FOR THE 

11 

2.4739E20 

1.2S50E18 

2.5208E18 


EXACT (16TH ORDER) 

AND APPROXIMATE (3RD 

12 

1.2973E21 

1.0446E19 

1.3685E19 


AND 9TH ORDER) 

TRANSFER FUNCTICTIS 

13 

4.2578E21 

5.1525E19 

4.6358E19 






14 

7.9983E21 

1.4079E20 

9.0249E19 


Poles of Poles of 

Poles of 

15 

7.4304E21 

1.7841E20 

8.7170E19 


16th order 3rd 

order 

9rh order 

16 

2.4119E21 

7.4230E19 

2.9263E19 


transfer 

Routh 

Routh 


functions 


approxinant approximant 


1 

TABLE II. IMPULSE RESPONSE ENERGY RATIOS FOR 2 

INCREASING APPROXIMANT ORDER 3 

(Fuel Flow Input) ^ 

5 

^ Routh Nc (p - 1, Tx (p - 2, 6 


approximant 

q ■ D- 

q ■ D- 

7 

order 

impulse 

Inipul.. 

8 

(k) 

response 

response 

9 


energy ratio 

energy ratio 

10 




11 

1 

0.41 

0.007 

12 

2 

.72 

.03 

13 

3 

.89 

.07 

14 

4 

.94 

.14 

15 

5 

.941 

.23 

16 

6 

.95 

.36 


7 

.96 

.51 


8 

.98 

.67 


9 

.99 

.81 


10 

.997 

.91 


11 

.998 

.97 


12 

.999 

.99 


13 

.9992 

.998 


14 

.9998 

.999 


15 

.9999 

.9999 


16 

1.0000 

1.0000 



-0.648 

-1.91 

-2.62 

-6.7111.31J 

-17.814.80J 

-18.2 

-21.611.55J 
-38.7 
-47.1 
-50.6 
-59.2 K 
-175. 

-577. 


-0.635 

-0. 98910. 580J 


-0.648 

-1.91 

-2.62 

-6. 5011. 16 J . 
-7.6413.29J 
-6.8218. 73J 
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